% This script is made to plot the occam velocity over
% occam predicted signals

% Sep 25, 2009 !

% Left=0;
% Right=360;
% Up = 90;
% Down=-90;
clear
Left=190;
Right=240;
Up = 60;
Down=45;

% Left=360-190;
% Right=360-160;
% Up = 60;
% Down=45;



occam = load('c:\manoj\ocean\occam_01_000\occam_01_000.txt');
L = (90-occam(:,1)) >= Down & (90-occam(:,1)) <= Up & occam(:,2) >= Left & occam(:,2) <= Right;
occam(~L,:) = [];
nlat = (Right-Left)/0.25;
nlong = (Up-Down)/0.25;
vel = load('C:\Manoj\ocean\OCCAM29853000\OCCAM29853000.txt');
L = (90-vel(:,1)) >= Down & (90-vel(:,1)) <= Up & vel(:,2) >= Left & vel(:,2) <= Right;
vel(~L,:) = [];
load C:\Manoj\m\figures\COLORTABLE.mat;

[Br] = reshape(occam(:,3),[nlat,nlong]);
[Bt] = reshape(occam(:,4),[nlat,nlong]);
[Bp] = reshape(occam(:,5),[nlat,nlong]);
[LatMatrix] = reshape(90-occam(:,1),[nlat,nlong]);
[LonMatrix] = reshape(occam(:,2),[nlat,nlong]);


figure(1);
set(gca,'FontSize',16)
worldmap([Down Up],[Left Right]);
geoshow(flipud(Br'),[4,LatMatrix(1,1),LonMatrix(1,1)],'DisplayType', 'texturemap');
h1=quiverm(90-vel(:,1),vel(:,2),vel(:,4),vel(:,3),3);
set(h1(2),'LineWidth',2);
set(h1(1),'LineWidth',1);
colormap(map);
caxis([-3,3]);
h=colorbar('horiz');
set(h,'FontSize',16)
geoshow('landareas.shp','FaceColor','none'); %
title('Br at sea level');
figure(2);
set(gca,'FontSize',16)
worldmap([Down Up],[Left Right]);
geoshow(flipud(Bt'),[4,LatMatrix(1,1),LonMatrix(1,1)],'DisplayType', 'texturemap');
quiverm(90-vel(:,1),vel(:,2),vel(:,4),vel(:,3),3);
colormap(map);
caxis([-3,3]);
colorbar('horiz');
geoshow('landareas.shp','FaceColor','none'); %
title('B_\theta at sea level');

figure(3);
set(gca,'FontSize',16)
worldmap([Down Up],[Left Right]);
geoshow(flipud(Bp'),[4,LatMatrix(1,1),LonMatrix(1,1)],'DisplayType', 'texturemap');
quiverm(90-vel(:,1),vel(:,2),vel(:,4),vel(:,3),3);
colormap(map);
caxis([-3,3]);
colorbar('horiz');
geoshow('landareas.shp','FaceColor','none'); %
title('B_\phi at sea level');



%to plot 1/12 deg OCCAM from Hemant
vel = load('c:\manoj\projects\global_hawk\200101_UV.xyz');
Left=160;
Right=210;
Up = 55;
Down=52.5;
L = (vel(:,2)) >= Down & (vel(:,2)) <= Up & vel(:,1) >= Left & vel(:,1) <= Right;
vel(~L,:) = [];
figure(1);
set(gca,'FontSize',16)
worldmap([Down Up],[Left Right]);
%quiverm(vel(1:10:end,2),vel(1:10:end,1),vel(1:10:end,4),vel(1:10:end,3),3);
quiverm(vel(:,2),vel(:,1),-vel(:,4),vel(:,3),3);
geoshow('landareas.shp','FaceColor','none'); %
title('Br at sea level');


